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ABSTRACT 

We numerically investigate the hydrodynamics of accretion disk reversal and relate 
our findings to the observed spin-rate changes in the accreting X-ray pulsar GX 1+4. In 
this system, which accretes from a slow wind, the accretion disk contains two dynami- 
cally distinct regions. In the inner part viscous forces are dominant and disk evolution 
occurs on a viscous timescale. In the outer part dynamical mixing of material with 
opposite angular momentum is more important, and the externally imposed angular 
momentum reversal timescale governs the flow. In this outer region the disk is split 
into concentric rings of material with opposite senses of rotation that do not mix com- 
pletely but instead remain distinct, with a clear gap between them. We thus predict 
that torque reversals resulting from accretion disk reversals will be accompanied by 
minima in accretion luminosity. 

Subject headings: accretion, accretion disks — binaries: general — hydrodynamics - 
methods: numerical — pulsars: individual (GX 1+4) — stars: rotation 
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1. Introduction 

The availability of continuous monitoring data 
for several X-ray pulsars as obtained with the 
BATSE All Sky Monitor on the Compton Gamma 
Ray Observatory over the last five years (Bild- 
sten et al. 1997, Nelson et al. 1997) has led to 
renewed interest in the spin evolution of accret- 
ing neutron stars with a strong magnetic field. 
From earlier pointed observations of X-ray pul- 
sars with other satellites it was already known 
that the spin rates of these neutron stars show 
changes on many timescales, ranging from unre- 
solved fluctuations (shorter than 1 day) to sys- 
tematic trends lasting over decades (see, e.g. the 
compilation by Nagase 1989). 

X-ray pulsars are usually divided into two 
groups, those which accrete from a Roche-lobe 
filling companion by way of an accretion disk, 
and those that accrete from the stellar wind of 
a companion star, in which case a disk may or 
may not be present. The second group tends to 
show more erratic spin rate changes on longer 
timescales, whereas the disk accreting pulsars 
generally show long term systematic trends. On 
the shortest timescales however the spin rate 
changes in the two groups appear to be compa- 
rable (de Kool & Anzer 1993, Baykal & Ogelman 
1993). 

Several models have been proposed to explain 
these spin rate fluctuations, or alternatively the 
torques that cause them. In the standard disk ac- 
cretion model (Ghosh & Lamb 1979) the pulsars 
are spinning at an equilibrium period, at which 
the spin-up torque due to the angular momen- 
tum of the accreting material (often referred to 
as the material torque) and the interaction of the 
magnetic field of the pulsar with the disk inside 
the co-rotation radius is balanced by a braking 
(spin-down) torque due to the interaction with 
the more slowly rotating outer parts of the accre- 
tion disk. The equilibrium period depends on the 
accretion rate, and the pulsar is expected to spin- 
up or down as the accretion rate increases or de- 



creases. This model appears consistent with the 
observed positions of X-ray pulsars in the pulse 
period - X-ray luminosity (equivalent to accre- 
tion rate) diagram. It also explains why the long 
term spin rate changes are considerably smaller 
than what would be expected from the material 
torque alone. 

However, the standard disk accretion model 
does not explain the spin rate fluctuations ob- 
served on very short timescales that are super- 
posed on the long timescale trends in most disk 
accreting sources. These fluctuations must be 
caused by torques that are as strong as the full 
material torque, much larger than the net torques 
expected to result from a small disequilibrium 
between spin-up and spin-down torques. 

For the wind-accreting pulsars, a different 
mechanism can be invoked to explain the positive 
and negative spin rate changes. Early numerical 
studies of Bondi-Hoyle accretion flow (e.g. the 
two dimensional simulations by Matsuda et al. 
1987, Fryxell & Taam 1988, 1989) demonstrated 
that it may be possible to form temporary ac- 
cretion disks with alternating senses of rotation 
in a wind accreting system. More recent high 
resolution three dimensional numerical investiga- 
tions of the so-called wind "flip-flop" instability 
may be found in Ruffert (1992, 1997) and the 
references therein. The timescale of the disk re- 
versals is found to be of the order of the sound 
crossing time of the accretion radius, which for 
the most common systems in which the neutron 
star accretes from a fast wind from an early type 
(O or Be) star is of the order of hours. This is 
consistent with the torque fluctuations observed 
in wind accreting sources, which also appear to 
be of the order of the material torque discussed 
above, but applied in random directions with the 
direction changing on a timescale shorter than 1 
day. 

The new continuous monitoring data obtained 
by BATSE (Bildsten et al. 1997, Nelson et al. 
1997, Chakrabarty et al. 1997) have yielded val- 
ues for the mass accretion rate (X-ray luminos- 
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ity) and the accretion torque (change in spin fre- 
quency) on a regular basis for a large number 
of X-ray pulsars. This data set has allowed a 
detailed comparison of the observed relation be- 
tween accretion rate and torque, and that pre- 
dicted by theoretical models. (It should always 
be kept in mind however that the fluxes measured 
by BATSE are only the pulsed fraction of the X- 
ray luminosity, and that it is not clear whether 
the ratio between pulsed and non-pulsed fluxes 
is constant). 

The observations are not consistent with the 
standard Ghosh & Lamb model, since this pre- 
dicts a clear correlation between spin-up and an 
increase in X-ray luminosity, whereas the obser- 
vations show a variety of behaviour, with spin- 
up or spin-down occurring at the same appar- 
ent luminosity. The cause of this type of be- 
haviour in disk accreting systems is not well un- 
derstood, but a possible explanation has recently 
been suggested by some of us (Li & Wickramas- 
inghe 1998). 

However, for wind accreting systems the possi- 
bility of accretion from temporary accretion disks 
with alternating senses of rotation is still ten- 
able. In analogy, Nelson et al. (1997) made the 
ad hoc suggestion that many observational fea- 
tures of some systems that are normally thought 
to contain disks (GX 1+4, 4U 1626-67) would 
be explained if they were accreting from disks 
with alternately prograde and retrograde senses 
of rotation. Previously, Makishima et al. (1988), 
Dotani et al. (1989) and Greenhill et al. (1989) 
had also sought to explain the rapid spin-down of 
GX 1+4 in terms of accretion from a retrograde 
disk. 

If the secondary star is feeding the accre- 
tion disk via Roche lobe overflow, as is almost 
certainly the case in 4U 1626-67, it is hard to 
conceive how a retrograde disk could ever come 
about. However, in the case of GX 1+4, the sug- 
gestion is not unreasonable. This X-ray pulsar is 
unique in the sense that it is accreting from a red 
giant or AGB star wind (Chakrabarty & Roche 



1997), and is in a very wide orbit. Estimating the 
timescale of disk reversal for accretion from such 
a wind, one obtains a timescale of the order of 
years, and the disk would form at a large radius 
(~ 10 13 cm) so that the inner part of the accre- 
tion flow is expected to be like a normal accretion 
disk. A timescale of years corresponds well with 
the timescale on which the accretion behaviour in 
GX 1+4 is observed to change, with a negative 
correlation between accretion rate and spin-up 
in some phases while the disk would be retro- 
grade, and a positive one at other times when it 
is prograde (Chakrabarty et al. 1997). Thus, this 
system is ideally suited to study the possibility 
of forming retrograde disks, since the timescale 
for disk reversal would be much longer than that 
of the torque fluctuations on a timescale of one 
day or less that are common in all types of X-ray 
pulsars. In the systems that accrete from a fast 
wind, the two timescales are comparable, and the 
effects will be difficult to separate. 

Relatively little modelling of the hydrodynam- 
ics of disk reversal has been performed, mak- 
ing it difficult to assess whether this model is 
in agreement with the observations. The previ- 
ous numerical studies of Bondi-Hoyle accretion 
flow (e.g. Ruffert 1997, Matsuda et al. 1987) did 
not resolve the inner part of the accretion disk, 
and did not yield much insight in such details as 
the interaction of an existing disk with material 
coming in with opposite angular momentum, the 
timescales determining the evolution of such sys- 
tems, and the relation between accretion torque 
and luminosity expected in such a model. 

In this paper we use Smoothed Particle Hy- 
drodynamics (SPH) simulations to start an in- 
vestigation into these processes, and compare our 
results with the observation of GX 1+4. 

2. Numerical Results 
2.1. Numerical Method 

The aim of the calculations described here is 
to investigate the response of an accretion disk 
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to a change in the specific angular momentum of 
the material it is being fed. In particular we are 
interested in the case were the specific angular 
momentum of the material added at the disk's 
outer boundary is reversed. 

The simulations were performed using an SPH 
code that has been specifically adapted for disk 
problems (Murray 1996). Two and three dimen- 
sional applications of the code can be found in 
Murray (1998) and Murray & Armitage (1998). 
To control the computational expense, motion 
perpendicular to the disk plane was neglected, 
and the calculations were run in two dimensions. 
We took full advantage of the adaptive nature of 
SPH and allowed the particle smoothing length 
(and hence the resolution) to vary in space and 
time. 

Viscous dissipation is included in the code us- 
ing a linear term in the SPH equations that is 
based upon the linear artificial viscosity term de- 
scribed in Monaghan (1992). Detailed descrip- 
tions and tests of the term's application to ac- 
cretion disks are to be found in Murray (1996, 
1998). In the continuum limit this dissipation 
term introduces a viscous force per unit mass 

a v = ^(V.(cSS)+V(cSV-v)), (1) 

where S is the deformation tensor, S is the col- 
umn density, h is the smoothing length, c is the 
sound speed, and £ is the SPH linear artificial vis- 
cosity parameter. In this work we have simply set 
C = 1- In the simulations described below, both 
the density and the sound speed vary on similar 
spatial scales to the velocity. Hence we cannot 
(as was previously done in e.g. Murray 1998) 
approximate the dissipation in equation |l] as the 
sum of a shear and a bulk viscosity. 

2.2. Simulation Details 

We simulated the gas flow in the annular re- 
gion r m < r < r m - y r; n j is the radius at which 
gas is added to the simulation. Ideally we should 
like r; n to be of the order of the magnetospheric 



radius r m (a radius characterising the truncation 
of the viscous accretion disc), and r in j to be of 
the order of 0.1 times the accretion radius, 

r acc = GM/vl, 

where G is the gravitational constant, M is the 
mass of the neutron star, and v w is the wind ve- 
locity. Indeed for sources accreting from a fast 
wind the disk formation radius is not much larger 
than r m and it is computationally feasible to sim- 
ulate the entire region r m < r < 0.1r acc . 

However, for the system GX 1+4, the accre- 
tion radius is several orders of magnitude larger 
than r m , and we can not resolve the full range 
in radius over which the disk is likely to extend. 
We have worked around this problem in the fol- 
lowing way. As discussed in the next section, we 
expect the disk to be divided in two parts, an in- 
ner one in which viscous forces dominate, and an 
outer one in which dynamical mixing of material 
with opposite angular momentum is more impor- 
tant. We therefore consider two cases, one with 
a very low viscosity which should apply to the 
outer regions of the disk, and another with much 
higher viscosity which models the parts of the 
disk where the transition from dynamical mixing 
to viscous transport occurs. Although this still 
does not cover the range in radius all the way 
down to the magnetospheric radius, we can ex- 
pect that inside the transition region a standard 
accretion disk is present which adjusts instanta- 
neously to its outer boundary condition. 

The gravitational potential was simply taken 
to be that of the accreting star, thereby neglect- 
ing tidal effects. This is likely a very good ap- 
proximation, since the hydrodynamical simula- 
tions of the flip-flop instability show that disk 
formation generally OCClirS Sit T ~ 0.1 Tacc ; 

and 

that r acc is always smaller than (or at most com- 
parable to) the Roche lobe radius. 

We chose a cylindrically symmetric mass in- 
flow condition at the outer boundary. Single par- 
ticles were added to the calculation at ten points, 
evenly spaced around a circle of radius cen- 
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tred on the neutron star. Each particle was set 
moving in the azimuthal direction with specific 
angular momentum j, corresponding to that of a 
Kepler orbit of radius r c i rc . Particles were added 
in this fashion at regular time intervals At. 

We did not attempt to model the details of 
the disk-magnetosphere interaction, and simply 
used a free inner boundary at r = r; n . At the 
end of each time step, any particle at a radius 
r < r m was assumed to have been accreted by the 
central object. The mass and angular momentum 
of the particle was recorded, and the particle was 
removed from the calculation. For computational 
reasons, particles at very large radii r > r ou t were 
removed in a similar fashion. However, in this 
calculation the mass (and angular momentum) 
lost at large radii was too small (less than 1% 
of the total angular momentum) to significantly 
effect the evolution of the disk. 

The temperature profile of a steady state disk 
was imposed upon the flow. Thus, the sound 
speed 

c = co(-^-r 3/8 (2) 

' circ 

where cq is the sound speed at the circularisation 
radius. We found that details of the disk's tem- 
perature profile did not significantly affect the re- 
sults described below. In other words, pressure 
effects were not important. 

In the simulation we have scaled quantities so 
that the mass of the accreting neutron star, M, 
the radius of mass addition, r; n j, and the Kep- 
lerian angular velocity at that radius, fiinj, are 
unity. The values of the various parameters used 
in this calculation are given in Table @ (in terms 
of our scaled units). 

2.3. Evolution with non-negligible viscos- 
ity 

We obtained the initial condition for this cal- 
culation by running the code for 200 Orj with 
mass added in an anti-clockwise direction as de- 
scribed above. At the end of this initialisation 
(time t = 0) the disk contained 16694 particles. 



The calculation proper was started by altering 
the mass addition routine so that newly added 
material flowed in the clockwise direction (with 
the same magnitude j as before). The simulation 
was then run for a further 200 fir. . 

The evolution of the disk's radial mass and 
angular momentum profiles from the commence- 
ment of clockwise mass addition are shown in 
Figs |l] and ^ respectively. The rates of mass 
and angular momentum accretion onto the cen- 
tral object are shown as functions of time in fig- 
ures H and For approximately the first 70 
fi^j of the calculation, the injection streams di- 
rectly impacted the disk, compressing it inwards. 
Fig |l] clearly shows the radial mass distribution 
becoming more strongly peaked over time, with 
the radius of maximum mass moving rapidly in- 
wards. Correspondingly the mass accretion onto 
the central object increased to a maximum at 
time t ~ 70 . This initial disk compression 
occurred on a timescale shorter than the outer 
disk's viscous timescale. The collision of the in- 
jection streams with the disk enforced a very ef- 
ficient exchange of angular momentum between 
the two flows. The specific angular momentum 
of the gas at the outer edge of the disk was thus 
sharply reduced. This material then moved to 
smaller radii, with radial pressure gradients and 
internal viscous torques too small to prevent the 
disk shrinking. 

Eventually, the disk was compressed to a size 
where the injection streams could no longer im- 
pact directly upon it. Consequently the disk 
shrank no further and the mass accretion dropped 
sharply away. 

It is easy to show that a ballistic particle, 
launched in the azimuthal direction at radius r; n j 
with the specific angular momentum of a Kepler 
orbit of radius r c ; rc , will reach periastron at ra- 
dius 



For our calculation, r p = n n j/3. Indeed, Fig [j] 
shows that the disk is compressed to this radius 
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at around t ~ 70 firj. The minimum compres- 
sion timescale is simply obtained by ignoring in- 
ternal stresses and calculating the angular mo- 
mentum that must be added to reduce the outer 
disk (plus newly added gas) to a circular orbit of 
radius r p . For this calculation At m ; n ~ 30 fij^j . 
In fact by t = 60 f2j^j all particles in the calcu- 
lation remaining at radii r > r; n j/3 already had 
their specific angular momentum reduced below 
that of Kepler orbit with radius ri n j / 3 . 

Meanwhile, the accretion streams were no longer 
hindered by the original disk, and were free to 
form a second disk. Initially a narrow ring 
formed at r = r c i rc . Both the outer and inner 
ring started to spread viscously, and the inner 
ring began to slowly accrete through r- m . Once 
the initially almost empty gap between the two 
rings filled via viscous diffusion, the very strong 
shear between the two disks allowed angular mo- 
mentum to be more efficiently extracted from the 
inner disk, and accelerated its demise. Hence the 
final drop in mass accretion was sudden rather 
than a gradual tailing away. 

The angular momentum accretion rate in this 
simulation was simply proportional to the mass 
accretion rate as is to be expected if gas at r- m is 
still in Keplerian motion. Not until the inner disk 
is entirely consumed does the accretion torque on 
the central object change sign. 

2.4. Evolution with negligible viscosity 

If viscous interaction in the disk is not impor- 
tant, the evolution of the disk must be different 
from the one discussed above, in which the be- 
haviour during the second half of the time inter- 
val modelled was dominated by viscous effects. 
We therefore performed further simulations with 
a much lower viscosity, obtained by reducing the 
sound speed by a factor of 5. The maximum 
smoothing length, /i max = 0.01 r; n j. As a result 
(see equation |l|) the viscous forces are approxi- 
mately a factor ten smaller than in the previous 
simulation. 



To save computation time, we used insight 
gained from the previous calculation and be- 
gan with two, well separated, oppositely rotat- 
ing rings with identical Gaussian density profiles 
centred at r c ; rc and 0.5r c j rc . New mass was then 
injected at the outer radius with a sense of ro- 
tation opposite to that of the outer ring, to sim- 
ulate a second reversal of rotation at the outer 
boundary. Mass resolution was improved by in- 
creasing the rate of particle injection by a factor 
5 over the first simulation. 

In fact we completed three separate simula- 
tions for this section in order to verify that the re- 
gions of interest were being adequately resolved, 
and that are results were not sensitive to the nu- 
merics. The first calculation began with 16826 
particles in the two rings combined. As in the 
previous section, particles were added at regular 
intervals from ten equally spaced points around a 
unit circle, with the points being fixed in the in- 
ertial frame. A second calculation was completed 
with identical system parameters but four times 
the mass resolution (the initial set up contained 
57680 particles). The rate of mass addition was 
identical to the first simulation. A third calcula- 
tion was done in which the mass resolution was 
identical to the first calculation but the mode of 
mass addition was varied. Particles were again 
added at ten equally spaced points around the 
unit circle, but each set of ten was randomly off- 
set in azimuth. This was done to confirm that 
the instability described below was real and not 
forced by the periodic nature of the stream im- 
pact points on the outer ring. We only identified 
minor differences in the results of the three sim- 
ulations. 

The evolution of the mass distribution for the 
third calculation is shown in figure As in the 
simulation described in the previous section, the 
outer ring moved in on the mixing timescale, but 
in this case it encountered the innermost ring. 
For this low-viscosity case, we found that the 
strong shear between the rings gave rise to a 
non-axisymmetric instability (see figure ||) that 
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deformed both rings severely. This instability 
mixed material with opposite angular momen- 
tum very effectively, so that the inner ring also 
fell in. The severe deformation was seen in 
all three simulations completed for this section. 
Once the outer ring reached the radius r p and was 
not forced in any longer, an empty gap between 
the rings formed, the rings rapidly became circu- 
lar again, and a new ring began forming at the 
circularisation radius. In figure [?] we have plot- 
ted particle neighbour number against radius for 
the third calculation at time t = 3 (by which 
stage the rings have largely recircularised and the 
gap between them has formed). At this time only 
260 of the 16864 particles at radii r < 0.3 r; n j 
have fewer than ten neighbours. This, together 
with the fact that similar results were obtained 
in the higher resolution simulation, indicates that 
the calculations were not under-resolved. 

We conclude from these calculations that if 
the externally imposed mass reversal timescale is 
significantly shorter than the viscous timescale at 
the circularisation radius, a number of concentric 
rings with alternating senses of rotation could be 
present inside the circularisation radius, down to 
the radius where the viscous timescale for the 
formation of an accretion disk or the interaction 
between two peaks becomes shorter than the re- 
versal timescale. 

3. Discussion 
3.1. Timescales 

From our results it is clear that there are two 
processes governing the inward transport of mat- 
ter through the accretion disk. One is the usual 
viscosity, and the other one is mixing of mate- 
rial with different (or opposite) specific angular 
momentum. A first estimate of the relative im- 
portance of these two processes in a physical sys- 
tem can be found by comparing the timescales 
on which they can change the disk structure. 

The viscous timescale for a standard disk is 
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in which ao.i is the Shakura-Sunyaev viscosity 
parameter in units of 0.1, v\q the radius in units 
of 10 10 cm, and h is the thickness of the disk. 
For a system like GX 1+4 this timescale is of 
the order of 7000 years if we assume that the 
disk forms at a radius of about 0.1 

^acc rN " / 
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cm. Note that counter-rotating rings will interact 
viscously on a timescale shorter than that given 
by the expression above, and can be estimated 
by replacing h/r by h/l, where I is the separation 
between the two peaks 

The typical timescale on which the inflow con- 
ditions at the outer boundary conditions can 
change due to the flip-flop instability is 



GM 
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(5) 



in which v w i^ is the wind speeds in units of 15 
km/s, and c S) io the sound speed in the wind in 
units of 10 km/s. The timescale needed by the 
dynamical processes associated with the change 
in boundary conditions to cause the disk struc- 
ture to change significantly must be the same 
as tff, since a significant change requires mix- 
ing material already present in the disk from the 
previous episode with a similar amount of mate- 
rial with different angular momentum. 

Again for the parameters of GX 1+4, we ex- 
pect a steady viscous disk to form when t v < tff, 
implying r < 5 x 10 10 cm if a and h/r have their 
typical values. For the pulsars accreting from a 
fast wind this limit corresponds to r < 10 8 cm, 
which is smaller than or comparable to the mag- 
netospheric radius so that in these cases viscous 
forces must play a negligible role in the accretion 
flow even if a disk forms. 

These considerations led us to examine two 
cases in the numerical simulations. In the first, 
the viscosity was so large that it governed the 
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evolution near the inner boundary of the compu- 
tational domain, which modelled the transition 
from the mixing dominated to the viscosity dom- 
inated regime at small radii. The second simula- 
tion addressed what happens in the outer parts of 
the disk where viscosity is unimportant. In both 
cases we found that it is always the externally 
imposed timescale that determined the timescale 
of changes at the inner boundary. If the viscos- 
ity was low, leading to a viscous timescale much 
longer than the external reversal time, dynam- 
ical effects took over to transport the material 
inward. This appears consistent with the long- 
term observed behaviour of GX 1+4. 

3.2. Mass and angular momentum accre- 
tion rates 

The mass and angular momentum rates shown 
in figures ^ and ^| vary relatively smoothly. This 
is to be expected as long as the disk is so large 
that its inner part is a steady viscous disk, which 
will smear out any more abrupt changes over a 
viscous timescale, which in turn is comparable 
to the externally imposed timescale as discussed 
above. These profiles should be fairly universal 
as long as viscosity is important close to the mag- 
netosphere, since the non-viscous evolution sim- 
ply pushes in a sequence of rings, with the inner 
one disappearing by forming a disk and being ac- 
creted as in our first simulation. 

The simulations indicate that a reversal of the 
material accretion torque should occur during a 
period of low luminosity. There may be some 
evidence for this in the data for GX 1+4, espe- 
cially the non-detection of the source when its 
long-time spin rate change reversed in 1984, or 
possibly the low luminosities just before and af- 
ter the 1994-1995 spin-up episode observed by 
BATSE (Chakrabarty et al. 1997). Although 
the timescales of these two observed features dif- 
fer by a factor of 10, it is not inappropriate to ap- 
ply the same model to them, since the timescale 
quoted in equation |^ for the change of the outer 
boundary condition is only typical, and shorter 



and longer timescales have been seen in calcula- 
tions of the flip-flop model. An alternative ex- 
planation put forward by Cui (1997) invokes a 
"magnetic propeller" mechanism to explain the 
low accretion episode of GX 1+4. 

We do not intend that this model explain all 
details of the observations of torque and luminos- 
ity in GX 1+4, which sometimes show the same 
torque for very different luminosities even when 
the disk is supposedly rotating in the same di- 
rection as evidenced by the short-term torque- 
luminosity correlation. Clearly an additional 
model for the interaction between disk and mag- 
netosphere and the associated braking torques is 
required. 

4. Concluding Remarks 

We are of course fully aware of the extreme 
simplifications that were used in our simulations. 
In more realistic studies, 3-dimensional effects 
such as subsequent bursts of material arriving 
at the outer boundary not being coplanar or not 
having exactly opposite angular momentum, or 
vertical motion due to the strong heating that 
must occur in the strong shearing regions should 
be included. We plan to address some of these in 
future work. 

We expect, however, that two main effects will 
remain standing even after more detail has been 
included. The first is that changes at the magne- 
tospheric radius occur on a similar timescale as 
changes at the outer boundary, and the second 
is that rings with opposite rotation do not tend 
to mix completely but rather create a gap be- 
tween themselves, leading to the prediction that 
torque reversals should occur during minima in 
accretion luminosity. 

The authors would like to thank Dayal Wick- 
ramasinghe and Phil Armitage for comments and 
discussions. 
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Fig. 1. — Mass as a function of radius for the high 
viscosity calculation (section 2.2), at the times 
shown. The angular velocity of mass added at 
the outer boundary is reversed at t = (topmost 
curve). Subsequent curves are offset 0.5 units in 
mass. Every second curve is dashed for clarity. 
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Fig. 2. — Angular momentum as a function of 
radius for the high viscosity calculation. Every 
second curve is dashed for clarity. 
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Fig. 3. — Mass accretion rate function of 
time for the high viscosity calculation. 
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Fig. 4. — Angular momentum accretion rate as a 
function of time for the high viscosity calculation. 
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Fig. 5. — Evolution of the radial mass profile 
for the third calculation with negligible viscosity 
(section 2.3) in which mass injection was "disor- 
dered" (i.e. the azimuth of each set of ten par- 
ticles was set randomly). The calculation begins 
(i = 0) with two counter-rotating Keplerian an- 
nuli with identical Gaussian density profiles cen- 
tred at r = 0.25 and r = 0.5r out . 
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Fig. 6. — Density maps of the "disordered in- 
jection" negligible viscosity calculation shown in 
figure |5[ In order to show as much detail in each 
frame, the maps have been separately scaled. 
The maximum density and time of each snap- 
shot are shown in the top and bottom right re- 
spectively of each panel. 
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Fig. 7. — Neighbour number plotted against ra- 
dius for each particle at t = 3 in the negligi- 
ble viscosity simulation with "disordered" mass 
injection. 
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Table 1 
Parameter values 



Parameter c h max r in r circ r out At 

Units r; n j f2j n j ?*inj Hiy ^inj ^*inj ^inj 

Value 0.1 ' 0.02 0.05 0.5 1.2 0.10 
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